In a book "Tsunami and Nonlinear Waves": 
Numerical Verification of the Hasselmann 
equation. 

Alexander 0. Korotkevich 1 , Andrei N. Pushkarev 2,3 , Don Resio 4 , and 
Vladimir E. Zakharov 5 ' 2,1 ' 3 

1 Landau Institute for Theoretical Physics RAS, 2, Kosygin Str., Moscow, 119334, 
Russian Federation kao@landau.ac.ru 

Lebedev Physical Institute RAS, 53, Leninsky Prosp., GSP-1 Moscow, 119991, 
Russian Federation 

3 Waves and Solitons LLC, 918 W. Windsong Dr., Phoenix, AZ 85045, USA 
andrei@cox . net 

4 Coastal and Hydraulics Laboratory, U.S. Army Engineer Research and 
Development Center, Halls Ferry Rd., Vicksburg, MS 39180, USA 

5 Department of Mathematics, University of Arizona, Tucson, AZ 85721, USA 
zakharov@math . arizona . edu 

Summary. The purpose of this article is numerical verification of the thory of weak 
turbulence. We performed numerical simulation of an ensemble of nonlinearly inter- 
acting free gravity waves (swell) by two different methods: solution of primordial 
dynamical equations describing potential flow of the ideal fluid with a free surface 
and, solution of the kinetic Hasselmann equation, describing the wave ensemble in 
the framework of the theory of weak turbulence. Comparison of the results demon- 
strates pretty good applicability of the weak turbulent approach. In both cases we 
observed effects predicted by this theory: frequency downshift, angular spreading as 
well as formation of Zakharov-Filonenko spectrum I u ~ lo~ 4 . To achieve quantita- 
tive coincidence of the results obtained by different methods one has to accomplish 
the Hasselmann kinetic equation by an empirical dissipation term Sdiss modeling the 
coherent effects of white-capping. Adding of the standard dissipation terms used in 
the industrial wave predicting model (WAM) leads to significant improvement but 
not resolve the discrepancy completely, leaving the question about optimal choice 
of S dl3S open. 

Numerical modeling of swell evolution in the framework of the dynamical equa- 
tions is affected by the side effect of resonances sparsity taking place due to finite 
size of the modeling domain. We mostly overcame this effect using fine integration 
grid of 512 x 4096 modes. The initial spectrum peak was located at the wave num- 
ber k = 300. Similar conditions can be hardly realized in the laboratory wave tanks. 
One of the results of our article consists in the fact that physical processes in finite 
size laboratory wave tanks and in the ocean are quite different, and the results of 
such laboratory experiments can be applied to modeling of the ocean phenomena 
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with extra care. We also present the estimate on the minimum size of the laboratory 
installation, allowing to model open ocean surface wave dynamics. 

1 Introduction. 

The theory of weak turbulence is designed for statistical description of weakly- 
nonlinear wave ensembles in dispersive media. The main tool of weak turbu- 
lence is kinetic equation for squared wave amplitudes, or a system of such 
equations. Since the discovery of the kinetic equation for bosons by Nord- 
heim [1] (see also paper by Peierls [2]) in the context of solid state physics, 
this quantum-mechanical tool was applied to wide variety of classical prob- 
lems, including wave turbulence in hydrodynamics, plasmas, liquid helium, 
nonlinear optics, etc. (see monograph by Zakharov, Falkovich and L'vov [3]). 
Such kinetic equations have rich families of exact solutions describing weak- 
turbulent Kolmogorov spectra. Also, kinetic equations for waves have self- 
similar solutions describing temporal or spatial evolution of weak - turbulent 
spectra. 

However, the most remarkable example of weak turbulence is wind-driven 
sea. The kinetic equation describing statistically the gravity waves on the 
surface of ideal liquid was derived by Hasselmann [4] . Since this time the Has- 
selmann equation is widely used in physical oceanography as foundation for 
development of wave-prediction models: WAM, SWAN and WAVEWATCH 
- it is quite special case between other applications of the theory of weak 
turbulence due to the strength of industrial impact. 

In spite of tremendous popularity of the Hasselmann equation, its valid- 
ity and applicability for description of real wind-driven sea has never been 
completely proven. It was criticized by many respected authors, not only in 
the context of oceanography. There are at least two reasons why the weak- 
turbulent theory could fail, or at least be incomplete. 

The first reason is presence of the coherent structures. The weak-turbulent 
theory describes only weakly-nonlinear resonant processes. Such processes are 
spatially extended, they provide weak phase and amplitude correlation on the 
distances significantly exceeding the wave length. However, nonlinearity also 
causes another phenomena, much stronger localized in space. Such phenom- 
ena - solitons, quasi-solitons and wave collapses are strongly nonlinear and 
cannot be described by the kinetic equations. Meanwhile, they could com- 
pete with weakly-nonlinear resonant processes and make comparable or even 
dominating contribution in the energy, momentum and wave-action balance. 
For gravity waves on the fluid surface the most important coherent structures 
are white-cappings (or wave-breakings), responsible for essential dissipation 
of wave energy. 

The second reason limiting the applicability of the weak-turbulent theory 
is finite size of any real physical system. The kinetic equations are derived 
only for infinite media, where the wave vector runs continuous D-dimcnsional 
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Fourier space. Situation is different for the wave systems with boundaries, e.g. 
boxes with periodical or reflective boundary conditions. The Fourier space of 
such systems is infinite lattice of discrete eigen-modes. If the spacing of the 
lattice is not small enough, or the level of Fourier modes is not big enough, 
the whole physics of nonlinear interaction becomes completely different from 
the continuous case. 

For these two reasons verification of the weak turbulent theory is an urgent 
problem, important for the whole physics of nonlinear waves. The verification 
can be done by direct numerical simulation of the primitive dynamical equa- 
tions describing wave turbulence in nonlinear medium. 

So far,the numerical experimentalists tried to check some predictions of 
the weak-turbulent theory, such as weak-turbulent Kolmogorov spectra. For 
the gravity wave turbulence the most important is Zakharov-Filonenko spec- 
trum ~ lo~ 4 [5]. At the moment, this spectrum was observed in numerous 
numerical experiments [6]- [19]. 

The attempts of verification of weak turbulent theory through numerical 
simulation of primordial dynamical equations has been started with numerical 
simulation of 2D optical turbulence [20], which demonstrated, in particular, 
co-existence of weak - turbulent and coherent events. 

Numerical simulation of 2D turbulence of capillary waves was done in [6] , 
[7], and [8]. The main results of the simulation consisted in observation of 
classical regime of weak turbulence with spectrum F w ~ w~ 19 / 4 , and discov- 
ery of non-classical regime of "frozen turbulence" , characterized by absence 
of energy transfer from low to high wave- numbers. The classical regime of 
turbulence was observed on the grid of 256 x 256 points at relatively high 
levels of excitation, while the "frozen" regime was realized at lower levels of 
excitation, or more coarse grids. The effect of "frozen" turbulence is explained 
by sparsity of 3- wave resonance, both exact and approximate. The classical 
regime of turbulence becomes possible due to nonlinear shift of the linear fre- 
quencies caused by enhanced level of excitation. Conclusion has been made 
that in the reality the turbulence of waves in limited systems is practically 
always the mixture of classical and "frozen" regimes. 

In fact, the "frozen" turbulence is close to KAM regime, when the dy- 
namics of turbulence is close to the behavior of intcgrable system [8] . 

The first attempt to perform modeling of the system of nonlinear waves 
(swell on the surface of deep ocean), solving simultaneously kinetic equa- 
tion and primordial dynamic equations, has been done in the article [15]. 
The results of this simulation again confirmed ubiquity of the weak-turbulent 
Zakharov-Filonenko asymptotic uj~ and shown existence of the inverse cas- 
cade, but presented essentially different scenario of the spectral peak evolu- 
tion. Detailed analysis shown, that even on the grids as fine as 256 x 2048 
modes, the energy transport is realized mostly by the network of few se- 
lected modes - "oligarchs" - posed in the optimal resonant condition. This 
regime, transitional between weak turbulence and "frozen" turbulence, should 
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be typical for wave turbulence in the systems of medium size. It was called 
"mesoscopic turbulence". Similar type of turbulence was observed in [17], [18]. 

In this article we present the results of new numerical experiments on 
modeling of swell propagation in the framework of both dynamical and kinetic 
equations, using fine grid containing, corresponding to 512 x 4096 Fourier 
modes. We think that our results can be considered as first in the world 
literature direct verification of wave kinetic equation. 

One important point should be mentioned. In our experiments we observed 
not only weak turbulence, but also additional nonlinear dissipation of the wave 
energy, which could be identified as the dissipation due to white-capping. To 
reach agreement with dynamic experiments, we had to accomplish the ki- 
netic equation by a phenomenological dissipation term Sdiss- In this article 
we examined dissipation terms used in the industrial wave-prediction models 
WAM Cycle 3 and WAM cycle 4. Dissipation term used in WAM Cycle 3 
works fairly, while Sdiss used in WAM Cycle 4 certainly overestimate nonlin- 
ear dissipation. This fact means that for getting better agreement between 
dynamic and kinetic computations, we need to take into consideration more 
sophisticated dissipation term. 



2 Deterministic and statistic models. 



In the "dynamical" part of our experiment the fluid was described by two 
functions of horizontal variables x,y and time t: surface elevation r](x,y,t) 
and velocity potential on the surface tp(x,y,t). They satisfy the canonical 
equations [23] 

d>q_5_H_ d^___5H_ 

' " SiP ' dt ~ ~ Sr) ' ( ' 



at 



Hamiltonian H is presented by the first three terms in expansion on powers 
of nonlincarity Vrj 



H = H 

H = 1 - 
2 
1 



H 1 +H 2 + 



Hi = 



H 9 



n 



dxdy, 



(2) 



r](kip) k(r)(kip)) + r]V 2, ip dxdy. 



Here k is the linear integral operator k = \J— V 2 , defined in Fourier space as 

^ J |£# k e- ikr dk, |fc| = 



klp r 



(3) 



Using Hamiltonian (2) and equations (1) one can get the dynamical equations 
[6]: 
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fj = ktp- (V^VVO) - Hvkip} + 
+k(r]k[rikil>]) + ±V 2 [r) 2 kip} + 

fktfV 2 ^] +F- 1 [ lkVk ] 1 ( 4 ) 

i> = -m-h [(v^) 2 -(^) 2 - 

— {kip}k[rjkip] — [rikipY^ 2 ^ + F -1 ^^]. 

Here i* 1-1 corresponds to inverse Fourier transform. We introduced artificial 
dissipative terms F _1 [7/ c V'fc], corresponding to pseudo- viscous high frequency 
damping. 

It is important to stress that we added dissipation terms in both equations. 
In fact, equation for r\ is just kinematic boundary condition, and adding a 
smoothing term to this equation has no any physical sense. Nevertheless, 
adding of this term is necessary for stability of the numerical scheme. 

The model (l)-(4) was used in the numerical experiments [6] - [8], [12], 
[13], [15], [17], [18]. 

Introduction of the complex normal variables ak 



ak= vS %+i vi^ k ' (5) 

where uj k = \fgk, transforms equations (1) into 

da k . 6H 

~~dt = (6) 

To proceed with statistical description of the wave ensemble, first, one 
should perform the canonical transformation a k — » 6k j which excludes the 
cubical terms in the Hamiltonian. The details of this transformation can be 
found in the paper by Zakharov (1999) [24]. After the transformation the 
Hamiltonian takes the forms 

H = J Wk&k&k + \J Tkkik^kkkAAa X ( ? ) 

X(5k+k 1 -k 2 -k 3 dkidk 2 dk3. 

where T is a homogeneous function of the third order: 

T(ek, eki, ek 2 , ek 3 ) = e 3 T(k, k 1; k 2 , k 3 ). (8) 

Connection between ak and 6 k together with explicit expression for T^^^ 
can be found, for example, in [24]. 

Let us introduce the pair correlation function 

< a k a k < >=giV k <5(k-k'), (9) 



where N-^ is the spectral density of the wave function. This definition of the 
wave action is common in oceanography. 
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We also introduce the correlation function for transformed normal vari- 
ables 

< Mk< >= gn^S{k - k') (10) 

Functions rik and iVk can be expressed through each other in terms of cumber- 
some power series [24]. On deep water their relative difference is of the order 
of /j, 2 (/x is the characteristic steepness) and can be neglected (in most cases 
of swell evolution (or wave evolution) experimental results shows /i ~ 0.1). 
Spectrum nk satisfies Hasselmann (kinetic) equation [4] 

—q^ = S n i[n] + S diss + 27 fc n k , 

S n i[n] = ^g 2 J |7k,k 1 ,k 2 ,k 8 | 2 (n kl n k2 n k3 + ^ 

+n k n k2 nk 3 - n^n^n^ - n^n^n^) x 

xS (ujk + Wfcj - Lu k2 - u)k 3 ) x 

xS (k + ki - k 2 - k 3 ) dkidk 2 dk 3 . 

Here Sdiss is an empiric dissipative term, corresponding to white-capping. 
Stationary conservative kinetic equation 

Snl = (12) 

has the rich family of Kolmogorov-type [25] exact solutions. Among them is 
Zakharov-Filonenko spectrum [5] for the direct cascade of energy 

n k ~ 1, (13) 

and Zakharov-Zaslavsky [26], [27] spectra for the inverse cascade of wave ac- 
tion 



3 Deterministic Numerical Experiment. 
3.1 Problem Setup 

The dynamical equations (4) have been solved in the real-space domain 2tt x 2tt 
on the grid 512 x 4096 with the gravity acceleration set to g = 1. The solution 
has been performed by the spectral code, developed in [21] and previously 
used in [22], [12], [13], [15]. We have to stress that in the current computations 
the resolution in Y"-direction (long axis) is better than the resolution in X- 
direction by the factor of 8. 

This approach is reasonable if the swell is essentially anisotropic, almost 
one-dimensional. This assumption will be validated by the proper choice of the 
initial data for computation. As the initial condition, we used the Gaussian- 
shaped distribution in Fourier space (see Fig. 1): 
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(15) 



The initial phases of all harmonics were random. The average steepness of 
this initial condition was \i ~ 0.167. 

To realize similar experiment in the laboratory wave tank, one has to to 
generate the waves with wave-length 300 times less than the length of the 
tank. The width of the tank would not be less than 1/8 of its length. The 
minimal wave length of the gravitational wave in absence of capillary effects 
can be estimated as \ m i n ~ 3cm. The leading wavelength should be higher 
by the order of magnitude A ~ 30cm. 

In such big tank of 200 x 25 meters experimentators can observe the evolu- 
tion of the swell until approximately 700To - still less than in our experiments. 
In the tanks of smaller size, the effects of discreetness the Fourier space will 
be dominating, and experimentalists will observe either "frozen", or "meso- 
scopic" wave turbulence, qualitatively different from the wave turbulence in 
the ocean. 

To stabilize high-frequency numerical instability, the damping function has 
been chosen as 



The simulation was performed until t = 336, which is equivalent to 92676, 
where T is the period of the wave, corresponding to the maximum of the 
initial spectral distribution. 

3.2 Zakharov-Filonenko spectra 

Like in the previous papers [10], [12], [13] and [15], we observed fast formation 
of the spectral tail, described by Zakharov-Filonenko law for the direct cascade 
n>k ~ fc~ 4 [5] (see Fig. 2). In the agreement with [15], the spectral maximum 
slowly down-shifts to the large scales region, which corresponds to the inverse 
cascade [26], [27]. 

Also, the direct measurement of energy spectrum has been performed dur- 
ing the final stage of the simulation, when the spectral down shift was slow 
enough. This experiment can be interpreted as the ocean buoy record - the 
time series of the surface elevations has been recorded at one point of the 
surface during T^uoy — 300T . The Fourier transform of the autocorrelation 
function 




(16) 



k d = 1024,7 = 5.65 x 10~ 3 . 



T bu „y/2 




(17) 
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allows to detect the direct cascade spectrum tail proportional to uo 4 (see 
Fig. 3), well known from experimental observations [28], [29], [30]. 
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Fig. 3. Energy spectrum in a double logarithmic scale. The tail of distribution fits 
to asymptotics u>~ 4 . 

3.3 Is the weak-turbulent scenario realized? 

Presence of Kolmogorov asymptotics in spectral tails, however, is not enough 
to validate applicability of the weak-turbulent scenario for description of wave 
ensemble. We have also to be sure that statistical properties of this ensemble 
correspond to weak-turbulent theory assumptions. 

We should stress that in our experiments at thebeginning |ak| 2 is a smooth 
function of k. Only phases of individual waves are random. As shows nu- 
merical simulation, the initial condition (15) (see Fig.l) does not preserve its 
smoothness - it becomes rough almost immediately (see Fig. 4). The picture of 
this roughness is remarkably preserved in many details, even as the spectrum 
down-shifts as a whole. This roughness does not contradict the weak-turbulent 
theory. According to this theory, the wave ensemble is almost Gaussian, and 
both real and imaginary parts of each separate harmonics are not-correlated. 
However, according to the weak-turbulent theory, the spectra must become 
smooth after averaging over long enough time of more than 1 //i 2 periods. Ear- 
lier we observed such restoring of smoothness in the numerical experiments of 
the MMT model (see [45], [46], [47] and [48]). However, in the experiments 
discussed in the article, the roughness still persists and the averaging does not 
suppresses it completely. It can be explained by sparsity of the resonances. 

Resonant conditions are defined by the system of equations: 

Wfe + LU kl = LU k2 + UJ k:3 , 

k + ki = k 2 +k 3 , 

These resonant conditions define five-dimensional hyper-surface in six-dimensional 
space k, ki,k2. In any finite system, (18) turns into Diophantine equation. 




(18) 
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Fig. 4. Surface |a k | 2 at the moment of time t ~ 67T . 



Some solutions of this equation are known [31], [17]. In reality, however, the en- 
ergy transport is realized not by exact, but "approximate" resonances, posed 
in a layer near the resonant surface and defined by 

\u3k + wfcj - Uk 2 - ujk+k x -k 2 \ < r, (19) 

where r is a characteristic inverse time of nonlinear interaction. 

In the finite systems k, ki, k2 take values on the nodes of the discrete grid. 
The weak turbulent approach is valid, if the density of discrete approximate 
resonances inside the layer (19) is high enough. In our case the lattice constant 
is Ak = 1, and typical relative deviation from the resonance surface 

~ <Ak ^^^-^2x lCr 3 . (20) 
ui u u 600 v ; 

Inverse time of the interaction r can be estimated from our numerical exper- 
iments: wave amplitudes change essentially during 30 periods, and one can 
assume: F juj ~ 10~ 2 » — . It means that the condition for the applicabil- 
ity of weak turbulent theory is typically satisfied, but the "reserve" for their 
validity is rather modest. As a result, some particular harmonics, posed in 
certain "privileged" point of fc-plane could form a "network" of almost reso- 
nant quadruplets and realize significant part of energy transport. Amplitudes 
of these harmonics exceed the average level essentially. This effect was de- 
scribed in the article [15], where such "selected few" harmonics were called 
"oligarchs" . If "oligarchs" realize most part of the energy flux, the turbulence 
is "mesoscopic", not weak. 
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3.4 Statistics of the harmonics 

According to the weak-turbulent scenario, statistics of the a^t) in any given 
k should be close to Gaussian. It presumes that the PDF for the squared 
amplitudes is 

P(|a k | 2 ) ~ I e -W 2 /^ ; (21) 

here D =< \a^\ 2 > — mean square amplitude. To check equation (21) we 
need to find a way for calculation of D(k). If the ensemble is stationary in 
time, -D(k) could be found for any given k by averaging in time. In our case, 
the process is non-stationary, and we have a problem with determination of 
D(k). 

To resolve this problem, we used low-pass filtering instead of time averag- 
ing. The low-pass filter was chosen in the form 

/(n) = e-d"!/- 4 ) 3 , A = 0.25iVx/2, Nx = 4096. (22) 

This choice of the low-pass filter preserves the values of total energy, wave 
action and the total momentum within three percent accuracy, see Fig. 5. Then 




Fig. 5. Low-pass filtered surface |ak| 2 at t ~ 67Tb. 



it is possible to average the PDF over different areas in fc-space. The results 
for two different moments of time t ~ 70T and t ~ 933T are presented in 
Fig. 6 and Fig. 7. The thin line gives PDF after averaging over dissipation 
region harmonics, while bold line presents averaging over the non-dissipative 
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Fig. 6. Probability distribution function (PDF) for relative squared amplitudes 
\a k \ 2 / < \a k \ 2 >. t~67T . 

area |k| < kd = 1024. One can see that statistics in the last case is quite close 
to the Gaussian, while in the dissipation region it deviates from Gaussian. 
However, deviation from the Guassianity in the dissipation region doesn't 
create any problems, since the " dissipative" harmonics do not contain any 
essential amount of the total energy, wave action and momentum. 

One should remember, that the bold lines in the Fig. 6 and Fig. 7 are the 
results of averaging over a million of harmonics. Among them there is a pop- 
ulation of " selected few" , or " oligarchs" , with the amplitudes exceeding the 
average value by the factor of more than ten times. The "oligarchs" exist 
because our grid is still not fine enough. 

In our case "oligarchs" do exist, but their contribution in the total wave 
action is not more 4%. Ten leading "oligarchs" at the end of the experiment 
are presented in the Appendix A. 

3.5 Two-stage evolution of the swell 

Fig. 8-11 demonstrate time evolution of main characteristics of the wave field: 
wave action, energy, characteristic slope and mean frequency. 

Fig. 10 should be specially commented. Here and further we define the 
characteristic slope as follows 
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Fig. 7. Probability distribution function (PDF) for relative squared amplitudes 
\a k \ 2 / < \a k \ 2 >. i~ 925T . 



According to this definition of steepness for the classical Picrson-Moscowitz 
spectrum \i = 0.095. Our initial steepness fi ~ 0.167 exceeds this value essen- 
tially. 

Evolution of the spectrum can be conventionally separated in two phases. 
On the first stage we observe fast drop of wave action, slope and especially 
energy. Then the drop is stabilized, and we observe slow down-shift of mean 
frequency together with angular spreading. Level lines of smoothed spectra in 
the first and in the last moments of time are shown in Fig. 12- 13 

Presence of two stages can be explained by study of the PDFs for eleva- 
tion of the surface. In the initial moment of time PDF is Gaussian (Fig. 14). 
However, very soon intensive super-Gaussian tails appear (Fig. 15). Then they 
decrease slowly, and in the last moment of simulation, when characteristics 
of the sea are close to Peirson-Moscowitz, statistic is close to Gaussian again 
(Fig. 16). Moderate tails do exist and, what is interesting, the PDF is not 




(23) 



Following this definition for the Stokes wave of small amplitude 



r\ = acos(kx), 
\i = ak. 
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Fig. 8. Total wave action as a function of time for the artificial viscosity case. 




Fig. 9. Total wave energy as a function of time for the artificial viscosity case 



symmetric — elevations are more probable troughs. PDF for 7] y — longitu- 
dinal gradients in the first moments of time is Gaussian (Fig. 17). Then in a 
very short period of time strong non-Gaussian tails appear and reach their 
maximum at t ~ 14T (Fig. 18). Here T = 2n / 'V^o — period of initial leading 
wave. Since this moment the non-Gaussian tails decrease. In the last moment 
of simulation they are essentially reduced(Fig.l9). 

Fast growing of non-Gaussian tails can be explained by fast formation of 
coherent harmonics. Indeed, 14Tq ~ 2tt /(^om) i s a characteristic time of non- 
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Fig. 10. Average wave slope as a function of time for the artificial viscosity case. 
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Fig. 11. Mean wave frequency as a function of time for the artificial viscosity case. 



linear processes due to quadratic nonlinearity. Note that the fourth harmonic 
in our system is fast decaying, Hence we cannot see "real" white caps. 
Figures 20-22 present PDFs for gradients in the orthogonal direction. 
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Fig. 15. PDF for the surface elevation r\ at some middle moment of time, t ~ 70Tn. 
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Fig. 16. PDF for the surface elevation r\ at the final moment of time, t ~ 933Tb. 
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Fig. 19. PDF for (Vr?) H at the final moment of time, t ~ 933T . 
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Fig. 20. PDF for (Vr/) x at the initial moment of time, t = 0. 
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Figures 23,24 present snapshots of the surface in the initial and final mo- 
ments of simulation. Fig. 25 is a snapshot of the surface in the moment of 
maximal roughness T — 4.94 ~ 14T - 
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Fig. 23. Surface elevation at the initial moment of time, t = 0. 
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Fig. 24. Surface elevation at the final moment of time, t ~ 933Tb. 
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Fig. 25. Surface elevation at the moment of maximum roughness, t ~ 14Tb. Gra- 
dients are more conspicuous. 
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4 Statistical numerical experiment 

4.1 Numerical model for Hasselmann Equation 

Numerical integration of kinetic equation for gravity waves on deep water 
(Hasselmann equation) was the subject of considerable efforts for last three 
decades. The "ultimate goal" of the effort - creation of the operational wave 
model for wave forecast based on direct solution of the Hasselmann equa- 
tion - happened to be an extremely difficult computational problem due to 
mathematical complexity of the S n i term, which requires calculation of a 
three-dimensional integral at every advance in time. 

Historically, numerical methods of integration of kinetic equation for grav- 
ity waves exist in two "flavors" . 

The first one is associated with works of [32], [33], [34], [35], [36] and [37], 
and is based on transformation of 6- fold into 3- fold integrals using S- functions. 
Such transformation leads to appearance of integrable singularities, which 
creates additional difficulties in calculations of the S n i term. 

The second type of models has been developed in works of [38] and [39] , [40] 
and is currently known as Resio- Tracy model. It uses direct calculation of res- 
onant quadruplet contribution into S n i integral, based on the following prop- 
erty: given two fixed vectors k, k x , another two k 2 , k 3 are uniquely defined 
by the point "moving" along the resonant curve - locus. 

Numerical simulation in the current work was performed with the help 
of modified version of the second type algorithm. Calculations were made on 
the grid 71 x 36 points in the frequency- angle domain [oj,0] with exponential 
distribution of points in the frequency domain and uniform distribution of 
points in the angle direction. 

To date, Resio- Tracy model suffered rigorous testing and is currently used 
with high degree of trustworthiness: it was tested with respect to motion 
integrals conservation in the "clean" tests, wave action conservation in wave 
spectrum down-shift, realization of self - similar solution in "pure swell" and 
"wind forced" regimes (see [42], [41], [43]). 

Description of scaling procedure from dynamical equations to Hasselman 
kinetic equation variables is presented in Appedix B. 

4.2 Statistical model setup 

The numerical model used for solution of the Hasselmann equation has been 
supplied with the damping term in three different forms: 

1. Pseudo- viscous high frequency damping (16) used in dynamical equations; 

2. WAMl viscous term; 

3. WAM2 viscous term; 

Two last viscous terms referred as WAMl and WAM2 arc the "white- 
capping" terms, describing energy dissipation by surface waves due to white- 
capping, as used in SWAN and W AM wave forecasting models, see [44]: 
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7k = C ds ul ((I -S)+ deltat) I J- ] (24) 

where k and u> are wave number and frequency, tilde denotes mean value; 
Cds, 5 and p are tunable coefficients; S = ky/H is the overall steepness; 
Spm — (3.02 x 10 -3 ) 1 / 2 is the value of S for the Pierson-Moscowitz spectrum 
(note that the characteristic steepness fi = ^/2S). 

Values of tunable coefficients for WAM1 case (corresponding to WAM 
cycle 3 dissipation) are: 

C ds = 2.36 x 1(T 5 , (5 = 0, p = 4 (25) 

and for WAM2 case (corresponding to WAM cycle 4 dissipation) are: 

Cds = 4.10 x 10~ 5 , (5 = 0.5, p = A (26) 

In all three cases we used as initial data smoothed (filtered) spectra (Fig. 5) 
obtained in the dynamical run at the time T* = 3.65rain = 24.3 ~ 70Xb, which 
can be considered as a moment of the end of the fist " fast" stage of spectral 
evolution. 

The natural question stemming in this point, is why calculation of the dy- 
namical and Hasselmann model cannot be started from the initial conditions 
(15) simultaneously? 

There are good reasons for that: 

As it was mentioned before, the time evolution of the initial conditions (15) 
in presence of the damping term can be separated in two stages: relatively fast 
total energy drop in the beginning of the evolution and succeeding relatively 
slow total energy decrease as a function of time, see Fig. 9. We explain this 
phenomenon by existence of the effective channel of the energy dissipation due 
to strong nonlinear effects, which can be associated with the white-capping. 

We have started with relatively steep waves /j, ~ 0.167. As we see, at that 
steepness white-capping is the leading effect. This fact is confirmed by nu- 
merous field and laboratory experiments. From the mathematical view-point 
the white-capping is formation of coherent structures - strongly correlated 
multiple harmonics. The spectral peak is posed in our experiments initially 
at k ~ 300, while the edge of the damping area k d — 1024. Hence, only the 
second and the third harmonic can be developed, while hire harmonics are 
suppressed by the strong dissipation. Anyway, even formation of the second 
and the third harmonic is enough to create intensive non-Gaussian tail of the 
PDF for longitudinal gradients. This process is very fast. In the moment of 
time T = 14Tb we see fully developed tails. Relatively sharp gradients mimic 
formation of white caps. Certainly, the "pure" Hasselmann equation is not 
applicable on this early stage of spectral evolution, when energy intensively 
dissipates. 

As steepness decreases and spectral maximum of the swell down-shifts, 
an efficiency of such mechanism of energy absorption becomes less important 
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when the steepness value drops down to \l ~ O.lthe white-capping becomes 
a negligibly small effect. It happens at T ~ 28070 ■ We decided to start com- 
parison between deterministic and statistical modeling in some intermediate 
moment of time T ~ 70T . 

5 Comparison of deterministic and statistical 
experiments. 

5.1 Statistical experiment with pseudo-viscous damping term. 

First simulation has been performed with pseudo-viscous damping term, 
equivalent to (16). 

Fig. 8 -11 show total wave action, total energy, mean wave slope and mean 
wave frequency as the functions of time. 

Fig. 32 shows the time evolution of angle-averaged wave action spectra as 
the functions of frequency for dynamical and Hasselmann equations. 

Temporal behavior of angle-averaged spectrum is presented on Fig. 32. We 
see the down-shift of the spectral maximum both in dynamic and Hasselmann 
equations. The correspondence of the spectral maxima is not good at all. 

It is obvious that the influence of the artificial viscosity is not strong 
enough to reach the correspondence of two models. 

5.2 Statistical experiments with WAM1 damping term 

Fig. 33 - 36 show total wave action, total energy, mean wave slope and mean 
wave frequency as the functions of time. 

The temporal behavior of total wave action, energy and average wave slope 
is much better than in the artificial viscosity term, and for 50 mm duration 
of the experiment we observe decent correspondence between dynamical and 
Hasselmann equations. However for longer time the WAM1 model deviates 
from the exact calculations significantly. 

It is important to note that the curves of temporal behavior of the total 
wave action, energy and average wave slope diverge at the end of simulation 
time with different derivatives, and the correspondence cannot be expected to 
be that good outside of the simulation time interval. 

Fig. 37 shows the time evolution of the angle-averaged wave action spectra 
as the functions of frequency for dynamical and Hasselmann equations. As in 
the artificial viscosity case, we observe distinct down-shift of the spectral max- 
ima. Correspondence of the time evolution of the amplitudes of the spectral 
maxima is much better then in artificial viscosity case. 
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5.3 Statistical experiments with WAM2 damping term 

Fig. 38 - 41 shows the temporal evolution of the total wave action, total energy, 
mean wave slope and mean wave frequency, which are divergent in time. 

Fig. 42 show time evolution of angle-averaged wave action spectra as the 
functions of frequency for dynamical and Hasselmann equations. While as 
in the artificial viscosity and WAM1 cases we also observe distinct down- 
shift of the spectral maxima, the correspondence of the time evolution of the 
amplitudes of the spectral maxima is worse than in WAM1 case. 

Despite the fact that historically WAM2 appeared as an improvement of 
W AMI damping term, it does not improve the correspondence of two models, 
observed in WAM1 case, and is presumably too strong for description of the 
reality. 

6 Down-shift and angular spreading 

The major process of time-evolution of the swell is frequency down-shift. Dur- 
ing T — 933T the mean frequency has been decreased from uj — 2 to u\ = .6. 
On the last stage of the process, the mean frequency slowly decays as 

< u> >~ r - 067 ~ i- 1 / 15 (27) 

The Hasselmann equation has self-similar solution, describing the evolu- 
tion of the swell n(k,t) ~ t 4 / n F (jtttt) (sec [41], [43]). For this solution 

<w>~t- 1/n (28) 

The difference between (27) and (28) can be explained as follows. What 
we observed, is not a self-similar behavior. Indeed, a self-silmilarity presumes 
that the angular structure of the solution is constant in time. Meanwhile, we 
observed intensive angular spreading of the initially narrow in angle, almost 
one-dimensional wave spectrum. Level lines of the spectra after low-pass fil- 
tering, obtained in dynamical equations simulation, for two moments of time 
are presented on Fig. 26-27. Level lines of the spectra in the same moments of 
time, obtained by solution of the Hasselmann equation are presented on Fig. 
28-29. One can see good correspondance between results of both experiments. 
Comparison of time-evolution of the mean angular spreading calculated from 
action and energy spectra are presented on Fig. 30-31. 




Fig. 26. Level lines of the spectra at t = 67Tb. Dynamical equations. 




Fig. 27. Level lines of the spectra at t = 674Tb. Dynamical equations. 




Fig. 28. Level lines of the spectra at t — 67Tq. Hasselmann equation. 




Fig. 29. Level lines of the spectra at t = 674Tb. Hasselmann equation. 
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Calculation from wave action spectrum 
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Fig. 30. Comparison of time-evolution of the mean angular spreading 
(/ |0|n(k)dk) / (Jn(k)dk) calculated through wave action spectra. 
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Fig. 31. Comparison of time-evolution of the mean angular spreading 
(J |#|um(k)dk) / (Jcjn(k)dk) calculated through wave energy spectra. 
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One has to expect that the angular spreading will be arrested at later 
times, and the spectra will take a universal self-similar shape. 
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Fig. 32. Angle- averaged spectrum as a function of time for dynamical and Hassel- 
mann equations for artificial viscosity case. 
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Fig. 33. Total wave action as a function of time for W AMI case. 



7 Conclusion 

1. We started our experiment with characteristic steepness /i ~ 0.167. This 
is three times less than steepness of the Stokes wave of limiting amplitude, 
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Fig. 34. Total wave energy as a function of time for W AMI case 
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Fig. 35. Average wave slope as a function of time for W AMI case. 

but still it is a large steepness typical for young waves. For waves of such 
steepness white-capping effect could be essential. However, in our experiments 
we cannot observe such effects due to the strong pseudo-viscosity. Indeed, third 
harmonics of our initial leading wave is situated near the edge of damping 
area, while fourth and higher harmonics arc far in the damping area. This 
circumstanse provides an intensive energy dissipation, which is not described 
by the Hasselmann equation. 

Anyway, on the first stage of the process we observe intensive generation of 
coherent higher harmonics which reveal itself in tails of PDF for longitudinal 
gradients. If our damping region would be shifted further to higher wave 
numbers, we could observe sharp crests formation. 
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Fig. 36. Mean wave frequency as a function of time for W AMI case. 
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Fig. 37. Angle- averaged spectrum as a function of time for dynamical and Hassel- 
mann equations a function of time for W AMI case. 

2. We ended up with steepness [i ~ 0.09. This is close to mature waves, typ- 
ically observed in the ocean and described by Hasselmann equation pretty well. 
We observed characteristic effects predicted by the weak-turbulent theory — 
down-shift of mean frequency formation, Zakharov-Filonenko weak turbulent 
spectrum w -4 and strong angular spreading. Comparison of time-derivatives 
of the average quantities shows that for this steepness wave-breaking (white- 
capping) become not essential at fi ~ 0.09. 

In general, our experiments validate Hasselmann equation. However, it has 
to be accomplished by a proper dissipation term. 

3. The dissipative term used in the WAM1 model fairly describe damping 
due to white capping on the initial stage of evolution. It overestimate damping, 
however, for moderate steepness [i ~ 0.09 
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Fig. 38. Total wave action as a function of time for WAM2 case. 
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Fig. 39. Total wave energy as a function of time for WAM2 case 

The dissipative term, used in the WAM2 model is not good. It overesti- 
mates damping essentially. 

4. Presence of abnormally intensive harmonics, so called "oligarchs" show 
that, in spite of using a very fine grid, we did not eliminated effects of discrete- 
ness completely. More accurate verification of the Hasselmann equation should 
be made on the grid containing more than 10 7 modes. This is quite realistic 
task for modern supercomputers, and we hope to realize such an experiment. 

Another conclusion is more pessimistic. Our results show that it is very 
difficult to reproduce real ocean conditions in any laboratory wave tank. Even 
a tank of size 200 x 200 meters is not large enough to model ocean due to the 
presence of wave numbers grid discreteness. 



Numerical Verification of the Hasselmann equation. 



35 



10 20 30 40 50 

Time (min) 

Dashed line - Hasselmann, solid line - dynamical equations 



Fig. 40. Average wave slope as a function of time for WAM2 case. 
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Fig. 41. Mean wave frequency as a function of time for WAM2 case. 
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9 Appendix A. "Forbes list of 15 largest harmonics. 

Here one can find 15 largest harmonics at the end of calculations in the frame- 
work of dynamical equations. Average square of amplitudes in dissipation-less 
region was taken from smoothed spectrum, while all these harmonics exceed 
level |a k | 2 = 1.4 x 1CT 12 . 
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9.1 Appendix B. Prom Dynamical Equations to 
Hasselmann Equation. 

Standard setup for numerical simulation of the dynamical equations (4), im- 
plies 2tt x 2ir domain in real space and gravity acceleration g = 1. Usage of 
the domain size equal 2ir is convenient because in this case wave numbers are 
integers. 

In the contrary to dynamical equations, the kinetic equation (11) is for- 
mulated in terms of real physical variables and it is necessary to describe the 
transformation from the "dynamical" variables into to the "physical" ones. 

Eq.4 are invariant with respect to "stretching" transformation from "dy- 
namical" to "real" variables: 

T) T = ar]' r ,, k=^-k', r = ar', g = vg' , (29) 

/ a 

t = W —t 1 , L x = aL' x , L y ~ a L' y (30) 

where prime denotes variables corresponding to dynamical equations. 

In current simulation we used the stretching coefficient a = 800.00, which 
allows to reformulate the statement of the problem in terms of real physics: 
we considered 5026 m x 5026 m periodic boundary conditions domain of sta- 
tistically uniform ocean with the same resolution in both directions and char- 
acteristic wave length of the initial condition around 22 m. In occanographic 
terms, this statement corresponds to the "duration-limited experiment" . 
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